#!/usr/bin/python

import sys
import re
from Bio.Blast import NCBIXML

anno_tab_file = sys.argv[1]
atf = open(anno_tab_file, "rU")
anno = atf.readlines()

m1 = re.compile('hypothetical protein')
m2 = re.compile('(pfam\d+), (\w+),.*')

for line in anno:
    l = line.split('\t')
    feat = l[0]
    locus_tag = l[1]
    start = int(l[2])
    stop = int(l[3])
    cds_length = (stop - start + 1) / 3
    frame = l[4]
    gene = l[5]
    protein_id = l[6]
    ncbi_ec_num = l[7]
    cog = l[8]
    bound_moiety = l[9]
    desc = l[10].rstrip()
    #ko = l[11].rstrip()
    outstr = ""
    if feat == "CDS":
        if m1.match(desc):
            #outstr = locus_tag
            #print locus_tag
            pass
        else:
            pfam_xml_file = locus_tag + ".pfam.rpsblast.xml"
            pxf = open(pfam_xml_file)
            pfam = NCBIXML.parse(pxf)
            pfam_rec = pfam.next()
            align_list = []
            hsp_lengths = []
            hsp_expects = []
            pfam_hit = ""
            annotation = ""
            if len(pfam_rec.alignments) == 0:
                outstr = locus_tag
            else:
                for alignment in pfam_rec.alignments:
                    if m2.match(alignment.hit_def):
                        pattern = m2.match(alignment.hit_def)
                        pfam_id = pattern.group(1)
                        domain = pattern.group(2)
                        align_list.append(domain)
                        #outstr = locus_tag + "\t" + pattern.group(2)
                        #total_hsp_length = 0
                        if len(alignment.hsps) == 1:
                            for hsp in alignment.hsps:
                                if hsp.expect < 1e-5:
                                    #total_hsp_length += hsp.align_length
                                    if hsp.align_length < 0.5 * alignment.length:
                                        print locus_tag + "\t" +  str(cds_length) + "\t" +  pfam_id + "\t" + domain + "\t" + "(" + str(hsp.align_length) + " < " + str(alignment.length) + ")"
                                    #hsp_lengths.append(hsp.align_length)
                                    #hsp_expects.append(hsp.expect)
                            #if total_hsp_length < 0.5 * alignment.length:
                            #    print locus_tag, cds_length, pfam_id, domain, total_hsp_length, "<", alignment.length
                        else:
                            total_alignment_length = 0
                            for hsp in alignment.hsps:
                                total_alignment_length += hsp.align_length
                                #if hsp.expect < 1e-5:
                                #    total_alignment_length += hsp.align_length
                            if total_alignment_length < 0.5 * alignment.length:
                                print locus_tag + "\t" + str(cds_length) + "\t" + pfam_id + "\t" + domain + "\t" + "(" + str(total_alignment_length) + " < " + str(alignment.length) + ")"
                    else:
                        outstr = locus_tag
                #outstr = locus_tag + "\t" + pfam_hit + "\t" + annotation
            pxf.close()
    #print outstr
atf.close()
